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ABSTRACT 

We present two-dimensional fully-kinetic Particle-in-Cell simulations of de- 
caying electromagnetic fluctuations. The computational box is such that wave- 
lengths ranging from electron to ion gyroradii are resolved. The parameters used 
are realistic for the solar wind, and the ion to electron mass ratio is physical. 
The understanding of the dissipation of turbulent fluctuations at small scales 
is thought to be a crucial mechanism for solar wind acceleration and coronal 
heating. The computational results suggest that a power law cascade of mag- 
netic fluctuations could be sustained up to scales of the electron Larmor radius 
and smaller. We analyse the simulation results in the light of the Vlasov lin- 
ear theory, and we comment on the particle heating. The dispersion curves of 
lightly damped modes in this regime suggest that a linear mechanism could be 
responsible for the observed steepening of power spectra at electron scales, but 
a straightforward identification of turbulent fluctuations as an ensemble of linear 
modes is not possible. 

Subject headings: 



1. Introduction 

The solar wind is a collisionless medium which is known to be turbulent over many 
length and time scales. The first observati ons of an active turbulent cascade at large scales 



in the solar wind date back to the work of lColemam (119681 ). In analogy with hydrodynamic 



turbulence, the large scale fluctuations are expected to nonlinearly cascade to smaller scales 
and being ultimately dissipated in the so-called dissipation range. However, in a collisionless 
plasma ordinary viscosity cannot play a decisive role for the dissipation of energy, as it does 



^now at the Applied Mathematics and Plasma Physics Group, Theoretical Division, Los Alamos National 
Laboratory, Los Alamos, NM 87545 



- 2 - 



in a neutral fluid. Hence, it has been long thought that the onset of dissipation range in 
the solar wind must be associated with cyclotron or Landau damping, and therefore must 
lie within kinetic scales. 

The understanding of the mechanisms that convert the electromagnetic turbulent energy 
into kinetic heating is far from complete. This issue is very relevant for the development of 
a realistic scenario of solar wind acceleration and coronal heating. 

The non-hyd ro dynamic nature of the physics involved in the dissipation process is now widely 
appreciated ^Smith et a/.l 120061 ) . and the understanding of dissipation mechanisms at small 
scales represents a very active field of research. 

Observati ons indicate that kine tic plasma physics at both proton and electron scales are 
involved ^Matthaeus et a/.ll2008[ ). Moreover, it is known that the Kolmogorov-like inertial 
range power spectrum terminates at ~ 0.3 Hz, where the spectrum of magnetic fluctuations 
steep ens to form a power la w with an average index close to —3 and varying between —2 and 
—5. nLeamon et a/.lll998al Jbl. 119991 ). More recent observations ^Sahraoui et a/.ll2009l . 120101 ) 
have shown that the power law in the range of scales between ion and electron gyroradii 
is quite universal, and close to 2.8. They also suggested the possibility that a nonlinear 
turbulent cascade could proceed further above the frequency associated with electron gyro- 
motion, where the power spectrum supposedly undergoes another abrupt break and steepen- 
i ng. T his interpretation has been challenged by observations reported in \Alexandrova et al 
( 120091 ). that suggest an exponential roll-over of the power spectrum at electron scales, which 
is usually interpreted as the onset of linear dissipative mechanisms. 

In the last few years an intensive effort has been directed to understand the mechanism 
that regulates the damping of turbulent fluctuations at small scales through a number o f 
different approaches involving numerical simulations. In particular, \Howes et al\ (j2008bl ) 
have performed state of the art gyrokinetic simulations of kinetic Alfven waves (KAW) cas- 
ca de; Particle- I n-Cel l (PIC) simulations of d ecaying whistler modes have been performed 
by \Saito et al\ (120081 ) ; \Svidzmski et al\ ( l2009l ) have presented P IC simulations of fast mag- 
netosonic waves. Hyb rid simulations have been presented by iParashar et al. (I20091) and 



Markovskii et a/.l (j2010[ ). and a shell model for 3D Hall-MHD has been studied by 



Galtier and Buchlin 



l2007f ). 

Interestingly, many works in this area have made the assumption that there must be a cer- 
tain length scale where the amplitude of magnetic fluctuations becomes so small that it is 
justifiable to treat the turbulent dissipation as damping of hnear waves. Therefore the spec- 
trum of magnetic fl uctuations has been tentatively compared to linear theory predictions in a 
number of papers (iLi et a/.ll2001t \Stawicki et a/.ll200ll : \Howes et a/.ll2008al : \ Gary et a/.l 12008 



numoer oi papers ii.Lz et aiwi^^ii \btawicKi et at.ll/uuil : \nowes et at.M/uugat \ uary et ai.M/uu»i : 
Jiang et a/.ll2009t \Podesta et al. koioh . Here it is (sometimes implicitly) assumed that linear 



damping mechanisms can convert electromagnetic energy to kinetic energy with the same 
efficiency and rate in a turbulent system and in an homogeneous plasma in thermodynamical 



- 3 - 



equilibrium. This approach has been recently challenged in lCamporeale et all (120 lOf ). where 
it is argued that using the least damped modes of a normal mode linear analysis might be a 
misleading way to characterise the dissipation. 

In any case, if any linear mode is dominant in regulating the casca de in the dissipation ra nge, 
there is certainly no agreement on which mode that should be ^Gary and SmitmuOO^ . It 
has to be mentioned that the numerical approach mentioned above and pursued by many 
researchers, has so far focused on ion scales, not extending very deeply into electron scales, 
mainly for reasons of computational costs. 

The aim of this paper is to present fully-kinetic PIC simulations of decaying electromagnetic 
fluctuations at electron scales. In doing so we address the following open problems. It is 
still unknown whether a linear approximation for the damping rate of the fluctuations, that 
would effectively describe the waves as an ensemble of linear modes, is fully justifiable and 
at which scale that approximation becomes valid. Moreover, there is an open debate about 
the scale at which such fluctuations are expected to be completely dissipated, and what 
mechanism is responsible for the dissipation. As we will show, our results do not support a 
scenario where turbulent fluctuations can be described as an ensemble of the least damped 
linear waves, not even at very small scales, although the observed steepening of the power 
spectra might be associated with the steepening of the dispersion curves of linear modes. 
The feasibility of the simulations presented in this paper has been greatly enhanced by using 
a semi-implicit method both for the field solver and the particle mover, and we will discuss 
some computational issues that have emerged in this study, and that should be taken into 
account for future studies. 

The paper is organised as follows. Section 2 briefly describes the PIC code, the parameters 
used, the issues related to the initialization of the simulations, and compares this paper with 
previous works. The main results are described in Section 3. The emphasis will be on the 
observed nonlinear cascade, the identification of linear modes, and particle heating. The 
conclusions and discussion of future work are given in Section 4. 



2. Methodology 

We perform simulations of an ion-electron plasma in a two-dimensional periodic box 
in Cartesian geometry {x,y). The code used is a fully-kinetic, electromagnetic, parallel 
PIC code, called PARSEK2D. The main feature of the code is the use of an implicit moment 
method for advancing the fields in time and a predictor-corrector routine to a dvance the par- 



ticles. A thorough description of the algorithm and the code can be found in \Markidis et al 



(120091 ). The implicitness of the scheme allows a relaxation of the stability conditions typical 
of an explicit code. In particular the cell size Ax is not required to be comparable with 
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the Debye length, the timestep At can be larger than the inverse plasma frequency and, in 
general, the factor cAt/Ax (where c is the speed of light) can be larger than one. 
For the simulations presented in this paper, we have found that it is still convenient, for ac- 
curacy reasons, to resolve the inverse electron plasma frequency. The cell size Ax, however, 
is about 20 times the Debye length, and the Courant factor cAt/Ax is about 9. Therefore 
the total saving with respect to an explicit code is at least a factor of 80,000 (not taking 
into account the complexity of the algorithm). This large saving factor translates into the 
possibility of simulating a plasma in realistic conditions in a computational box which is 
fairly large. We compare in Figured] the range of wavevectors included in the computational 
box for the simulations presented in this paper with previous works. Table [1] compares other 
relevant parameters, such as the ratio of plasma to cyclotron frequency (this is undefined in 
hybrid simulations), and the ion to electron mass ratio. Notice that the use of an artificial 
mass ratio makes the values of the ion and electron gyroradii closer to each other, resulting 
in an extension of the range of wavevector studied. This of course results in a computa- 
tional advantage, but it is often difficult to estimate the effects of this approximation on 
the results, since it is very well known that the separation of ion and electron scales is an 
intrinsic and important characteristic of plasmas. What emerges from Figure [T] is that, as 
far as we know, this paper presents the first simulations that use realistic parameters in a 
fairly large c omputation a l dorn ain that extends to electron scales. We point out that also 
the work by \Saito et all (120081 ) was in principle able to extend to electron scales, but the 



number of particles used was too low in order to distinguish the results above kpi ~ 10 from 
the noise level. 

In this work, the plasma is in solar wind conditions. The electron plasma beta is equal to 0.5 
(/3e = 87inoTe/ Bq, Te is the electron temperature, no is the density, and Bq is the magnetic 
field intensity). Ions and electrons have initially equal temperature, and a Maxwellian dis- 
tribution function. The ratio of ion plasma frequency to ion gyrofrequency Upi/Qi is about 
1650, and the ion to electron mass ratio is physical (mj/me = 1836). The number of cells 
is 200 X 200, and the box size is large enough to accommodate waves with wavenumbers k 
ranging from kpi = 4.28 (equivalent to kpe = 0.1), to kpe = 10, where pi and pe are the ion 
and electron gyroradius, respectively. In all the simulations presented in this paper we have 
used 6400 particles per cell, for each species, and the simulations are typically run for a time 
T = 100r2~^ (fig is the electron gyrofrequency). 

The choice of initial conditions offers a wide degree of freedom. The general idea is to ini- 
tialize a spectrum of fiuctuations with properties somehow similar to the one observed in the 
solar wind, and follow their decay in time. Some authors have performed simulations where 
a spectrum of only one plasma linear mode was excited at small values of k (for instance, 
KAW or whistlers). In order to do so, one should compute the phase and amplitude relations 
between different quantities such as magnetic and electric field, velocity and density pertur- 
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bation predicted by linear theory for any single value of k, and impose such perturbations to 
the background equilibrium. In order to get a completely self-consistent initial perturbation, 
the particle distribution function should be perturbed exactly. This is very expensive, from 
a computational point of view, and a shortcut which is very often employed is to initialize 
the particles with a shifted Maxwellian distribution, such that only the current and charge 
densities are consistent with the linear theory predictions, but not all the higher-order mo- 
ments of the distribution function. The implicit assumption of such a method is that the 
PIC code will be able to quickly 'self-adjust' the initial inconsistency, but this approach 
remains, in our view, not very satisfactory. Moreover, the choice of the slope of the initial 
power spectrum constitutes yet another degree of freedom in the initialization. 
In this work, we have opted for a very general initial perturbation. In this way, we do not 
face the problem of choosing some particular modes to initialize, and we are also able to 
study which mode, if any, 'survives' through the nonlinear cascade. Accordingly, we initial- 
ize our simulations with a background magnetic field in the x direction, and we superimpose 
a spectrum of random phase fluctuations in the magnetic field only. We assume that such 
perturbation is general enough to perturb a large variety of waves in the plasma. We point 
out that, differing from the normal mode initialization described above, this choice is self- 
consistent both with Vlasov and Maxwell equations. The plasma will certainly respond 
instantaneously to magnetic field gradients by creating currents and generating electric fluc- 
tuations, and ultimately waves. The only requisite that the initial perturbation must satisfy 
is to be divergence-free. Hence, we initialize the magnetic perturbations as: 
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where Lx and Ly are the box length, kx = 2'K'm/Lx^ ky = 27rn/Ly, m and n are integers 
which usually range from —3 to 3, and (px^y, <Pxi <Pyi 4'z are random angles. Note that the 
choice of 6Bz is arbitrary and is made in such a way that its amplitude is of the order of the 
amplitudes of 6Bx and 6By. We have run several simulations with different values for the 
slopes a and /3 of the spectra and the results seem not to be significantly affected by those 
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parameters, For simplicity, we will show the results when a = /3 = 0. 

As a further check on the soundness of the results we have also initialized some simulations 
with a spectrum of waves with an Alfvenic velocity perturbation imposed, and the results 
are qualitatively similar, for the length and time scales investigated here. 

3. Results 

3.1. Turbulent cascade 

The first results that we show are for benchmark cases, where we have initialised the 
perturbation [SB /Bo ~ 1) with a single wavenumber and not with a spectrum of waves. Fig- 
ure [2] shows the developed spectrum of magnetic fluctuations 6B/Bq in wavenumber space, 
at the end of the simulation. The only mode initially perturbed was the mode m = 1, n = 3, 
which is indicated with an arrow. The result indicates that the harmonics of the initial mode 
get predominantly excited according to a 3-wave resonant interaction {k = ki + k2)- Clearly 
one cannot speak of turbulent cascade in this case, but the purpose of this simulation is to 
check whether the code is able to treat in a consistent and satisfactory manner injections 
of energy at small scales. An issue with PIC codes is indeed the unavoidable presence of 
noise at small scales due to the discreteness of computational particles. This noise could 
result in the spurious effect of an inverse cascade of energy from large wavenumbers which 
is completely numerically generated. One can see that this effect is not important in Figure 
[21 where the modes at large k that get excited are only the higher harmonics of the initial 
mode. The price to pay in order to avoid spurious effects at small scales is the use of an 
unusually large number of particles (6400 particles per cell). 

A confirmation of the fact that 3-wave resonant interactions are well described in the code 
is shown in Figure [21 In this case we have excited, along with the mode m = l,n = 3, its 
counter-propagating mode m = — l,n = 3. The modes initially excited are circled. The 
resulting distribution of energy satisfies again the relation k = ki + k2, giving rise to that 
checkerboard plot. 

As we have said, these results should only be considered as benchmark runs. However, from 
Figure [21 and [31 it is already clear that the treatment of small amplitudes fluctuations at 
large k as linear modes poses some problems when the system allows the presence of finite 
amplitude perturbations at small k. This is because nonlinear interactions that couple small 
and large k modes could always take place, resulting in a continuous injection of energy at 
small scales. Hence, in order to characterise the dissipation range solely through some linear 
damping mechanism, one needs to find out at what scale the nonlinear coupling becomes 
completely ineffective. In other words, the fact that the amplitude of fluctuations is small 
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at some scales does not justify alone the validity of linear theory, in the solar wind. Another 
requisite is that those scales should not take part in any nonlinear coupling with larger scales. 
We now show and comment on the results of a run where a fiat spectrum of modes has been 
excited, accordingly to Eqs. dlM]). The range of modes excited is (0, 3) in both wavenumbers 
m and n. As we said, we have tried different slopes for the initial spectrum, but we will 
comment on one run only, as the results were very similar. 

Figure m and [5] show the spectra of |(5i?p/|i?oP in parallel and perpendicular directions (i.e. 
integrated along the other direction), respectively. The red curve indicates the noise level, 
which is the level of fluctuations generated by particle noise, when no initial perturbation is 
imposed in the equilibrium. As it is clear, results above kpe ~ 7 are not reliable, because the 
curves approach the noise level. Moreover, in order to be far enough from the noise level, 
and allow the cascade to proceed as much as possible, we initialise the perturbation with a 
relatively large amplitude {6B/Bq ~ 3), at kp^ = 0.1. This is not completely consistent with 
solar wind observations, but it does not affect the purpose of this study, because the ampli- 
tude of fluctuations above kpe = 1 is small enough, in principle, to justify the comparison 
with linear theory which is presented in the next section. The spectra presented in Figures 
m and E] are fully converged. The cascade proceeds very quickly and reaches a stationary 
profile after about TQ^ ~ 30. Of course, in reality, there may be evolution on time scales 
longer, and length scales larger than those simulated. This could only be investigated with 
considerably larger and longer simulations. 

Two interesting features emerge from Figure H] and [51 First, the spectrum develops a well 
defined power law that steepens at about the inverse electron gyroradius. The plots show 
that the least-squared fit with a power law function yields exponents which are consistent 
with solar wind observations (in particular \6B/Bo\^ ~ /c~^'^ between the ion and electron 
inverse gyroradii, although this feature might be coincidental due to the 2D character of 
the simulations). Second, the nonlinear cascade proceeds far into the electrons scales, until 
it reaches the noise level. There is no sign of an exponential roll-over or the emergence of 
linear dissipation mechanisms, even when the level of SB/Bq is as low as 10"^, which usually 

considered small enough for applying linear theory. 

The nonlinear cascade proceeds as it is expected in magnetized plasma a Chen et a/.ll2010l ). It 
develops the well-known power anisotropy, which means that most of the energy resides, at 
the end of the simulation, in oblique or quasi-perpendicular modes. This is visible in Figure 
ini where we show the spectrum of magnetic perturbations as a contour plot in wavevector 
space. 
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3.2. Comparison with linear theory 

In previous works it has been argued that the dynamics of the dissipation at small 
scales can be described by linear mechanisms, since the amplitude of fluctuations become so 
small that non-linear effects can be neglected. It is also argued that a purely linear damping 
mechanism will result in an exponential roll-over of the spectrum of turbulent fluctuations. 
In this section, we attempt to interpret the simulation results at the light of Vlasov linear 
theory. We have used a Vlasov solver to identify linear modes in the range of wavevectors 
of interest for this study. We anticipate that, at these scales, it becomes very difficult 
to follow the numerical solution of the dispersion relation for one single mode, as many 
branches of the solutions overlap. We show in Figures [31S] the dispersion relations for whistler 
modes, KAW and Langmuir waves, for different angle of propagation {9 = 45°, 60°, 80°). 
We focus attention on oblique modes, since we know that the cascade tends to develop 
the aforementioned anisotropy. We point out that what we call Langmuir waves are the 
generalization of e l ectros tatic Langmuir modes in presence of magnetic field, studied in 



Willes and Cairnm ( l2000l ). They are lightly damped high-frequency oscillations, which are 
not usually taken in account in this regime. 

The same dispersion curv es for the damping rate only of KAW and whistler are plotted in 



log-log scale in Figure HI \Li et all (120011 ) have proposed an empirical fit for the damping 
rate 7 of KAW of the form: 

7 = -mi/t^^e-C^)', (4) 

where mi,m2,m3 are positive fitting parameters. They have also argued, on the basis of 
a simplified model, which treats the nonlinear cascade as an isotropic diffusion process in 
wavenumber space, that in order to recover a power-law for the spectrum of magnetic fluctu- 
ations, the damping rate should be itself a power-law as a function of k. What is interesting, 
in Figure ini is that at this scale the damping rate seems to follow precisely a power law in k. 
The form of Equation H] could still be valid, but the exponential factor becomes essentially 
very close to 1, as A; becomes very large. 

Another interesting feature of Figure which has not been previously commented on else- 
where, is that the damping rate of KAW and whistler undergo a smooth transition from a 
certain power law to another. In principle, one could argue that the dispersion curves of 
Figure [7] and |9] might account for the results of Figures H] and O i.e. for the cascade of the 
power spectrum and its steepening, since the properties of the linear waves, at this scale, 
seem to be favourable for the support of a power law solution for the fluctuations. However, 
if this is the case, we still need a firm characterisation of the fluctuations in terms of linear 
waves. 



Following \Gary and Smitm fj2009[ ) we now try to characterize the linear waves by their elec- 



tron compressibility 6nf./SB. We prefer to use the electron compressibility rather than the 
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magnetic compressibility, because the former gives a better signature of the difference be- 
tween KAW and whistler modes. Indeed, for KAWs the electron compressibility grows as a 
function of k, while it decreases for whistlers. In Figure [10] we show the values of the electron 
compressibility calculated from the simulation for different angle of propagation (diamonds 
for 6 = 45°, stars for 6 = 60°, and circles for 6 = 80°). The solid lines indicates the value 
predicted for KAWs, whistler and Langmuir modes. The simulations results clearly indi- 
cates a trend of increasing values for larger values of kpe, which suggests a certain degree 
of Alfvenicity of the fluctuations. However, none of the linear modes taken in consideration 
is in good match with the simulation result. The fluctuations are probably a mixture of 
different modes. Interestingly, the Langmuir mode seem to be the one whose electron com- 
pressibility is closest to the simulations results. This might be related to the high level of 
initial fluctuations. 



3.3. Particle heating 

The energy at the beginning of the simulations is partitioned equally between the mag- 
netic field and the particles, and electrons and ions have the same kinetic energy. The 
simulation conserves more than 95% of the total energy, and at the end of the run there is a 
5% increment of the total kinetic energy, at the expenses of magnetic energy. The increase 
in electric energy is negligible. The heating is directed mainly to electrons, such that at the 
end of the simulation they are 20% hotter than the protons. This is however not surprising, 
given the time and length scales of the simulations. 

Electrons are heated predominantly in the parallel direction (with respect to the magnetic 
field). The value of temperature anisotropy is not constant over the whole computational 
domain, but depends on the spatial position, and T\\/T_i_ ranges from 0.98 to 4.9. This sug- 
gest that the heating is connected with the formation of magnetic topological structures. 
An interesting feature that emerges by examining the electron distribution function is the 
presence of non-thermal features. However, the deviation from a bi-Maxwellian distribution 
is more accentuated when the particles and the mean magnetic field are considered in a small 
portion of the box. We show in Figure (ITT!) the electron distribution function in the (x, y) 
plane for four different nested boxes. Velocities are normalized to the speed of light. The 
solid line indicates the direction of the mean magnetic field within each box, and L is the 
length of the boxes, in ion gyroradius units. It turns out that as the magnetic field and the 
particles are averaged over larger and larger boxes, the distribution function tends to become 
more and more bi-Maxwellian and any non-thermal feature is lost. A similar result holds in 
the [x, z) plane, while the distributions in the direction orthogonal to the mean magnetic 
field remains roughly gyrotropic. This result has important implications for satellite obser- 
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vations. Indeed, it suggests that the time window over which the particles and the magnetic 
field are averaged affects the distribution functions and therefore the anisotropy measured. 
By assuming a steady wind flowing at 250 km/s, one can transform the box lengths into fre- 
quencies at which data are collected, and estimate what is the level of anisotropy measured 
at different frequencies. The upper bound of such anisotropies are shown, for electrons, in 
Figure (fT2l) . Given the small length of the computational box used in this study, the range 
of frequencies plotted in Figure (fT2|) is relatively high. However, future missions and instru- 
ments may be able to collect data at such frequencies, and therefore will be able to verify 
this prediction. 



4. Conclusion 

We have performed two-dimensional fully-kinetic PIC simulations of decaying magnetic 
fluctuations at electron scales. The scope of this work is to address some open questions 
regarding the dissipation of turbulent fluctuations at small scales, and the plausibility of 
treating such fluctuations in the framework of linear theory. 

We have shown the results of the flrst simulations that use realistic parameters in a fairly 
large computational domain that extends to electron scales. This represents an improvement 
in the context of the previous 2D simulations. 

The results of our simulations suggest that a turbulent cascade in the form of a power-law 
could proceed very far below the electron gyroradius. It is interesting that at such scales the 
damping rate for both kinetic Alfven waves and whistler modes are power laws as a function 
of k. In principle, this would suggest that purely linear mechanisms could be responsible 
for a power law spectrum of fluctuations, such as the one observed in satellite data, and 
obtained in the PIC simulations. However, we have shown that neither KAW or whistler 
have an electron compressibility which is in agreement with simulations results. 
Our interpretation is that many different modes should be taken in account at this scale, 
including Langmuir waves and possibly other more heavily damped modes. Unfortunately, 
the level of excitation of each mode depends on non-linear mechanisms active at larger scale. 
Therefore we doubt that the dissipation process can be understood and described solely 
within a linear description. However, the puzzling steepening of the power spectrum might 
be caused by changes of the linear properties of the waves, as the ones reported in Figure HI 
Concerning the heating of the particles, we have observed a preferential heating of electrons 
along the parallel direction (the same happens, in smaller measure, to protons, but the size 
of the box does not really allow to draw any informative conclusions on protons). It has 
also emerged that non-thermal features in the electron distribution function appear when 
the box over which quantities are averaged is not too large (i.e. smaller than 0.2-0.3 pi). 
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On the basis of this result we have made a prediction of the upper bound value of electron 
anisotropy related to the frequency at which data are collected, which can be verified by 
future measurements. 

A last comment on the computational aspects of this work is due. The results presented in 
this paper have been obtained essentially thanks to the use of a semi-implicit code. This has 
allowed to use a computational box which comprehends scales as little as kpe = 10 and as 
large as kpi ~ 4. However, in order to keep a low level of noise we had to use a large number 
of particles per cell. At present, this excludes the possibility of extending the box so to fully 
include ion scales, and constrains the geometry to two-dimensional. We used a 2D compu- 
tational spatial domain which contains the background magnetic field; the particle velocity 
space is 3D. This restricts the available /c— space for fluctuations to a (A;||,A;_l) plane, with 
the consequence, particularly, of restricting, but not entirely eliminating the interaction of 
modes with k strictly perpendicular to the magnetic field. However, the use of this particular 
choice of simulation plane retains fluctuations with a k\\ component and the possibility of 
particle parallel scattering and heating, and is thus a more realistic choice than having the 
background magnetic field strictly perpendicular to the simulation plane. For these reasons, 
we believe the 2D results we have presented are qualitatively robust, but, of course, they 
will only be confirmed by fully three-dimensional simulations. Of course, this is now at the 
edge of computational feasibility, but a possibility which is worth investigating in the future 
is to perform an implicit 5/— PIC simulation, so that the number of particles needed could 
be reduced. 
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Table 1: Comparison of different parameters with previous works. p.j is the ion gyroradius, 
uji/Vti is the ratio of ion plasma to cyclotron frequency (a typical value for the solar wind 
at 1 AU is around 4000), mj/mg is the ratio of ion to electron mass, Np is the number of 
particles per cell. 
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Fig. 1. — Comparison of the range of wavevectors studied in different works. The acronym 
list is reported in Table [H The vertical dashed lines indicate the inverse proton Larmor 
radius {kpi = 1), and the inverse electron Larmor radius {kpe = 1), assuming equal ion and 
electron temperatures. Notice that the latter is shifted in works that do not use a physical 
mass ratio. Also, notice that VIO has reported simulations with different temperature ratios 
(and therefore different locations of kpe = 1). 
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Fig. 2. — Spectrum of magnetic fluctuations 5B/Bq in wavenumber space, in logarithmic 
scale. The mode initially excited m = l,n = 3 is indicated with an arrow. 
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Fig. 3. — Spectrum of magnetic fluctuations 5B/Bq in wavenumber space, in logarithmic 
scale. The modes initially excited m = 1, n = 3 and m = —1, n = 3 are circled. 
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Fig. 4. — Spectrum of magnetic fluctuations |5i?p/|5oP in the parallel direction kpe\\. The 
noise level curve is in red. The power law best fits are superimposed. 




Fig. 5. — Spectrum of magnetic fluctuations |(5i?p/|i?oP in the perpendicular direction kpe±- 
The noise level curve is in red. The power law best fits are superimposed. 
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Fig. 6. — Contour plot of the spectrum of magnetic fluctuations 
perpendicular directions. The contours are for values of 10~^^, 10~^^, . . 
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.,10° 




Fig. 7. — Dispersion curves for KAW and whistler for different angles of propagation: 9 = 
45°, 60°, 80°. Frequencies are normalized on the electron cyclotron frequency fl^ (^^r is the 
real frequency and 7 is the damping rate). 
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Fig. 8. — Frequency (top) and damping rate (bottom) for the Langmuir mode at 6' = 80°, 
normalized on the electron plasma frequency Upe- 




Fig. 9. — Damping rates for KAW and whistler for different angles of propagation: 6 = 
45°, 60°, 80°, in log-log scale, normalized on the electron cyclotron frequency Qe- 



- 20 - 




Fig. 10. — (Color online) Electron compressibility |(5?7,e|/|(5i?|. The symbols are the results 
of simulations (circle for 9 = 80°, star for 6 = 60°, diamond for 6 = 80°). Solid lines are the 
predictions from linear Vlasov theory for KAW, whistler and Langmuir waves (in blue for 
^ = 80°, red for ^ = 60°, and black for ^ = 80°). 
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Fig. 11. — (Color online) Electron distribution functions in the {x,y) plane, collected in four 
nested boxes of increasing size L. The solid line shows the direction of the mean magnetic 
field within each box. Velocities are normalized to the speed of light. 
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Fig. 12. — Electron anisotropy T\\/T^ versus the frequency at which data are collected, in 
Hz. 



